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Abstract. 

A finite Larmor radius approximation is derived from the classical Vlasov equation, in the limit of large (and 
uniform) external magnetic field. We also provide an heuristic derivation of the electroneutrality equation 
in the finite Larmor radius setting. Existence and uniqueness of a solution is proven in the stationary frame 
for solutions depending only on the direction parallel to the magnetic field and factorizing in the velocity 
variables. 

Introduction. The ITER project is a challenge to the growing need of new sources of energy. It aims at 
producing energy by nuclear fusion. Nuclear reactions take place in a tokamak, where a high temperature 
plasma is confined. So far, confined plasmas are performed with relatively short energy confinement times 
due to microscale instabilities that generate turbulent transport |III99) . It is observed that the characteristic 
frequencies of these instabilities is several orders of magnitude smaller than the ion Larmor gyration frequency 
governed by the strong magnetic field. Studies of nuclear fusion in tokamaks are in full expansion, both 
experimentally and theoretically. Kinetic models are appropriate for studying the core of the plasma since 
the collisions have a very weak effect in these hot and low density plasmas. Physicists use gyrokinetic models 
and especially the finite Larmor radius approximation to model the core of the plasma [GIVT09 . Taking 
into account the fast Larmor gyration of the charged particles that characterizes magnetic confinement, these 
models allow one to average over that fast gyration and reduce the 6D kinetic problem to a 5D gyrokinetic 
one. That property is especially interesting for numerical simulation. In this paper the finite Larmor radius 
approximation is derived from the Vlasov equation, in the limit of large uniform magnetic field and with an 
external electric field. Because of the homogenization on the fast Larmor gyration, the limit equation (jl.lOp 
is written in 5D gyro-coordinates (x g ,v\\, \vj_\) defined in (jl.llj) . These coordinates are the position of the 
so-called particle guiding center, x g in the 3D space together with the parallel velocity v» and the amplitude 
of transverse velocity \v± \ that is proportional to the magnetic moment, an adiabatic invariant of the particle 
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motion in the strong magnetic field limit (given a constant magnetic field) . Its mathematical structure is a 
combination of the Vlasov equation in the direction parallel to the magnetic field and of the Euler equation 
in the perpendicular direction, where the original helds are replaced by the corresponding gyro-average fields. 

To close the system, physicists use the electroneutrality equation n e = Zn^ where n e stands for the 
electron density and rii the ion density, each ion having a charge Z . In the following Z = 1 will be considered 
with no loss of generality. It can be shown that the electroneutrality equation is in fact the the Poisson 
equation solved on scales that are significantly larger than the Debye length. Since the latter governs the 
Laplacian term of the Poisson equation, the electroneutrality equation is an appropriate approximation on 
scales of the order or larger than the Debye length. Moreover, taking into account the difference between the 
density of particles and that of guiding centers leads one to introduce a polarization correction due to the 
non uniform distribution of particles on gyro-circles. In the second part of this paper, the electroneutrality 
equation (|I.25I) is carefully written. The coupling of the 5D gyrokinetic Vlasov equation (|1.I0[) to the 
electroneutrality equation (|1.25[) is the model used for instance in the GYSELA code, a project that aims 
at modeling the turbulent transport in fusion plasmas [GSG + 06bl |GSG + Q6a . While in the GYSELA code 



there is the possibility to use a collision operator we concentrate here on the actual Vlasov equation with no 
collisions. 

A difficulty raised by the electroneutrality equation taken as such is the lack of an explicit regularization 
term for the electric potential. Consequently, the regularity of the latter is not sufficient to ensure a mathe- 
matical solution of the Vlasov transport equation. The analysis of the well-posedness of limit model based on 
the electroneutrality equation thus seems impossible today. For that reason we restrict ourselves to the II? 
Vlasov equation (without gyro-kinetic effects), coupled to the electroneutrality condition. That model is a 
particular case of the 3D model, provided the solutions do not depend on the direction perpendicular to the 
magnetic field. Even for that simple kinetic model, not much is known about solutions. The difficulty lies 
in the fact that the force term is proportional to the derivative of the density in the field direction. Indeed, 
while the Vlasov equation ensures that all L p -norms of the distribution function can be bounded, there is no 
control on the norms of its derivatives. In this paper, we prove the existence and uniqueness of a steady state 
solution in a slab geometry, therefore between two boundary conditions, provided one fulfills some conditions, 
in particular that they are no particles trapped between the two boundaries. 



1 Derivation of the finite Larmor radius 
approximation 

The ion distribution in low density plasmas submitted to a magnetic field is well described by the Vlasov 
equation. The latter is valid provided one can neglect the two-particle distribution function altogether 
[Nic83 so that the evolution of the standard distribution function is governed by the Liouville equation for 
the conservation of particles 

% + v ■ V x f + —(E(t, x) + v x B) • V„/ = , (1.1) 
at m% 

where / is the phase space density depending on time, position (in the domain D of the plasma) and velocity 
(in K 3 ) of the ions, and where Ze and are the ion charge and mass respectively. As stated in the 
introduction, we shall consider Z = 1 in the following with no loss of generality. A priori, the electric and 
magnetic fields are governed by the Maxwell equations, but a scale analysis allows one to approximate, and 
thus simplify, these laws. This will be made clear in the following. 

For a strong external magnetic field, the charged particles exhibit a fast rotation motion around the 
magnetic field lines. The frequency of that gyration, the Larmor frequency, is several orders of magnitude 
larger than the observed frequency range of turbulence. Furthermore, in the quasineutral limit this frequency 
is larger than all other frequencies of the particle dynamics, thus providing the means for an efficient scale 
separation. It is noteworthy that this is the basis of magnetic confinement that is presently realized in devices 
such as the tokamaks. In this framework, one separates two parts in the particle motion, on the one hand 
the slow motion of the center of the Larmor gyration, and, on the other hand, the fast gyro-motion. The 
phase space reduction achieved by only considering the slow motion, thus ignoring the fast gyro-motion, 
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leads one to the 5D gyrokinetic model. These are usually derived by physicists either by averaging the single 
particle dynamics ( [Cat 78] and [CTB81 ) or by a Lie-transform perturbative approach QHah.88] . jBH07j and 
[GLV08]). Mathematically, this is achieved by looking at the limit of equation (jl.ip when the modulus of 
the external magnetic field |-B ex t| tends to infinity. However, different models can be obtained, depending on 
the way the other control parameters vary when \B ext \ — > oo. For a magnetized plasma, two relevant scales 
characterize the limit regimes that have been discussed. These are the Debye length Ad and the thermal 
Larmor radius of the ions pth- The Debye length weights the Laplacian in the Poisson equation and thus 
defines the transition from the non-neutral description of the plasma required on the sub-Debye scales from 
the quasineutral plasma description for scales larger than the Debye length. The thermal Larmor radius is the 
radius of the fast Larmor gyration for ions with average speed Vth- The averaging over the small time scales 
governed by the ion Larmor frequency then tends to translate into an averaging over the ion Larmor scale. 
The finite Larmor radius terms are therefore introduced to take into account this rather weak cut-off effect. 
The radius of the electrons gyration is generally much smaller for comparable ion and electron temperatures 
and will be neglected. The two scales discussed above are defined by 



where the electron thermal energy T e is introduced rather than that of the ions since it is more appropriate 
to characterize the large electron mobility and therefore strong response to the electric field. Let us introduce 
the characteristic scale of the system L, for instance introduced by the boundary conditions. The large 
magnetic field limit then corresponds to the vanishing ion Larmor radius limit, namely the gyro-center 
approximation where pth/L — > 0. Within this framework, E. Grenier uses pseudo-differential calculus to 
prove }Gre97] the convergence of a 2D fluid model towards an incompressible fluid model when e goes 
to zero with pth/L = e and (Ad/L) 2 = e. For the same scaling and cold initial distributions, i.e. the 
approximation of a Dirac measure at velocity zero for the distribution function, Y. Brenier proved in [BreOOj 
with modulated energy technics the convergence of solutions of the Vlasov-Poisson system towards dissipative 
solutions (introduced by P.-L. Lions) of the Euler equation, for weel-prepared initial data. F. Golse and L. 
St-Raymond used the same technics in |GSR03] to prove the convergence of a Vlasov-Poisson model on the 
torus towards the 2D1/2 Euler equation. This is done for a vanishing s with pth/L — e and using another 
scaling for the Debye length, namely Xd/L = e. This property is restricted to a torus size such that there 
is no resonant oscillation. This particular case appears to be relevant for electrons in a region close to the 
tokamak boundary. 

Derivations have also been performed in the gyrocenter approximation, pth/ L±_ small but finite. Here Lj_ 
is a carcteristic length in the perpendicular direction, which is usually choosen smaller thant the one in the 
parallel direction in finite Larmor radius approximation. In [FSOOa] (see also |FS00b| for a short version) E. 
Frenod and E. Sonnendriicker studied the convergence of a linear Vlasov equation (external magnetic field), 
in this finite Larmor radius limit (large B but finite Larmor radius), using two scale convergence methods. 
M. Bostan studied the 2D strong magnetic field limit, using Hilbert expansion technics, in a setting of finite 
Debye length and Larmor radius, \d/{eL) = pth/{sL) — 1, where s is small but finite [Bos09] and obtained 
strong convergence results for regular initial conditions. Here, we further examine the linear case studied by 
Frenod and Sonnendriicker, and obtain a simpler and more complete derivation. The finite Larmor radius 
approximation is correct for fusion plasmas, including the core and most of the edge plasma. However, the 
ion Larmor radius is much larger than the Debye length throughout the plasma. An interesting derivation 
should be done in the framework of present gyrokinetic calculations, namely pthj{sV) = I, \d/{eL) — > 0, 
yet nothing has been done in this very demanding limit (see the beginning of If .21 for more details). 

The resolution of equation (|f .![) is known in the case where the fields E and B are external and C 1 or 
at least BV . In the BV-c&se, one can use the DiPerna Lions theory of transport equation p3L89 b] with its 
latest developments ( Amb04] ) . In the case where the electric potential is given by self-interaction with the 
usual Poisson law, and the magnetic field B is still considered as external, results of existence and uniqueness 
obtained for the Vlasov equation without magnetic field may be used with the appropriate modifications. 
We refer to |LP91j . |Hor93j . In the case where the self-induced magnetic field is not negligible any more, 
we refer to the work [DL89a of DiPerna and Lions about the Vlasov-Maxwell equation. In the case of an 
electric field given by an electroneutrality equation, nothing is known from the mathematical point of view. 




Pth = 



eB 



(1.2) 
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1.1 The finite Larmor radius gyro-kinetic approximation. Rigorous approach. 

For the sake of simplicity, we neglect the variation of the external magnetic field B ext . In this cylindrical 
approximation, the curvature of the magnetic field is not considered, and one neglects the exploration of the 
magnetic field variation by the particles during the fast cyclotron motion. With respect to the relevant scales 
this effect if of the order of the Larmor radius p t h divided by the characteristic radius of the magnetic field 
curvature, namely the tokamak major radius R. Although pth/R is a small parameter, this approximation 
is too strong with respect to many aspects of the tokamak physics, in particular since curvature effects are 
considered by physicists as the cause of a large class of micro-instabilities. The strong impact of such a small 
parameter can be traced back to a symmetry breaking of the system since the curvature governs the leading 
term that produces a charge separation. The present simplification must be regarded as a first step that is 
only valid when one addresses issues that do not lead one to symmetry breaking. Moreover, in order to avoid 
problems with boundary conditions, which are complex to handle, a periodic setting is used. In other words, 
we work on the torus T 3 in space, and R 3 in speed, with a constant magnetic field B = |£?|(0,0, 1). 
Next, we define dimensionless variables 

,/ t i x \\ i x x , v i E 
t = -, x = — , x x = — , V = — , E = — , 
t 11 L n L± v t h E 

with characteristic time t, parallel length Ln, perpendicular length L±, velocity Vth and electric field Eq. In 
this new set of dimensionless variables the Vlasov equation (II. lj) may be rewritten as: 

df . Tv th Tv th _ , reE reB ± 

7^77 + -r-d x i f+ —— ■ V x < / E ■ V„// ■ V // = , (1.3) 

of L|| II L± ± miV th m l 

where the subscript _L (resp. ||) stands for the projection on the perpendicular (resp. parallel) direction, 
and the superscript _L stands for the projection on the plane orthogonal to the field line and the rotation by 
— 7r/2 around the field line direction: v 1 - — (v2,— v\,Q) if v = (vi,V2,V3). The time scale r can conveniently 
be defined to reduce the number of control parameters in Eq. dl.3p by setting: 

rvth _ 1 

In a similar fashion, the normalization of the electric field can be used to define the third control parameter, 
hence: 

TeE _ 
miVth 

The two remaining control parameters are then L\\/L±_ for the third term in Eq. dl.3p and Ln/pth for the 
last term. It will be assumed here that both parameters L± and pth exhibit the same asymptotic behavior 
L± cx pth- Let us define the perpendicular scale as Lj_ = pth^TT / \h_\_Pth) ■ In this expression, k± is the typical 
wave vector of the cross-field fluctuations. For a sufficiently small Larmor radius, one can then assume that 
the turbulence follows the so-called gyro-Bohm scaling such that k±p t h is a constant LEHT02 . In particular 
it does not depend on the magnitude of the magnetic field. One then finds the following relationship between 
the two control parameters 

Ln L» 1 1 



Pth Lj_ k^Pth £ 

that can be grasped by the following relation whenever one drops the proportionality constant 

L \\ L \\ 1 
pth L± e 

In practice, one could also define the normalization scale as Lj_ = p t h so that the the two control parameters 
would be exactly the same. In this discussion, we have considered a transverse scale related to the fluctuation 
properties that allows one to characterize one of the terms in the Vlasov equation. However, other control 
parameters must be considered to account for the boundary conditions. The so-called p* parameter used in 
magnetic fusion is such a parameter since = pth/a where a is the plasma minor radius, the scale related to 
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the boundary conditions while pth is the fluctuation scale. Relating the parameter e to p* introduces aspect 
ratios that stem from the periodic boundary conditions. 

a 




The ratio a/Ln depends of the safety factor and the tokamak aspect ratio in the actual tokamak geometry but 
represents the ratio of the sizes of the domain in the radial and parallel direction with the present setting. 
The parameter e can also be expressed in terms of the slow and fast characteristic times of the particle 
motion. The slow time r introduced in the normalization of Eq. (jl.3[) . is the characteristic time to explore 
the tokamak geometry along the parallel direction and the fast time is due to the gyration motion, hence: 

1 

In this last step one thus finds that the ion Larmor frequency f2 = eB/rm is of order (re) -1 . In other 
words, one thus assumes that the ion Larmor frequency is much larger than the parallel connection frequency 
r . That is a valid assumption in all regions of a tokamak. Indeed, for ITER conditions, the ion Larmor 
frequency is of the order of 2 10 s Hz and the connection time r ranges from 10~ 3 s to 10 _4 s. These parameters 
mainly depend on the magnitude of the magnetic field and on the size of the device so that e characterizes 
a given fusion device, e ~ 10 -5 for ITER parameters. 

Dropping the primes, we obtain the following rescaled version of the Vlasov equation 

fit i 

-KT+v ■ d Xu f + E-V v f+ -(«! • V X J + v x ■ V vx f) = , (1.4) 



When e — > 0, the largest terms are those proportional to the factor 1/e. Retaining the latter, we obtain a 
transport equation associated to the following system of ODE in the plane transverse to the magnetic field, 

1.1, , , 

x = — v±, v = —v ± . (1-5) 

e e 

As B is homogeneous, the trajectories are circles of center x g — x± + and radius \v±\ covered with 
frequency e _1 . The global motion is the sum of this very quick motion of gyration and a slower and more 
complicated motion (with velocity of order one). In the limit of large \B\, particles are assumed to be evenly 
distributed on the gyro-circles, and their motion is the sum of a drift in the perpendicular direction, and a 
classical acceleration in the parallel direction. Heuristically, the electric drift ve may be obtained from the 
Newton law in normalized variables, 

v = E+-v x , (1.6) 

£ 

Given the assumption that the particles have a fast motion of gyration, one can integrate the previous 
equation over a period of gyration, at lowest order, ve is given by the averaged velocity such that: 

&{E) + - ve- so that v E = £ (E x ) . (1.7) 

In the latter expression one readily recognizes the usual form of the electric drift velocity ve = E x B/B 2 , 
however where the electric field is averaged over a period (E). One can show that this average of the field 
E corresponds to an average over a circle of radius \v±\ in the perpendicular plane. Provided the only 
dependence on the gyrophase stems from the particle motion, this average translates into a Bessel operator 
defined as: 

J° PL h ( x 9 ) = ^J h(x g +p L e^)d Vc , (1.8) 

where p L is the ion Larmor radius and e l¥>c = (cos (p c , sinyj c , 0). We also introduce the operator which 
stands for a position-velocity version of this average, 

Jp L 9(x g ,P L ,v l{ ) = —j^ g(x g + p L e tlfi %p L e l ^-i) +vn ) d <p c , (1.9) 
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where en = (0, 0, 1). 

Given the average on the gyrophase ip c , the motion is reduced to a 5D space. The gyration radius 
p L = \v±\ (given the chosen normalization) is related to the magnetic moment p, — mip 2 L /(2eB), a quantity 
that is an adiabatic invariant, i.e. a constant of motion in the large B limit. One can show that this invariant 
is the conjugate variable of the gyrophase. Here, as B is homogeneous, p L — |t>jj will remain constant in the 
limit \B ext \ — > +oo. We thus obtain a AD + ID model, i.e. a 5D model with no dynamics in the variable 
derived from the magnetic moment. This reduction of the phase space is very interesting for numerical 
simulations. The following theorem states this reduction precisely. 

Theorem 1.1 Let us assume that f® is uniformly bounded in L q , q > 1, that it weakly converges toward 
f G L q , and that E is a gradient and belongs to L\ x where p^ 1 + q^ 1 = 1. For each e > 0, let f £ be a 
solution of (|1.4[) with the initial condition f®. Then, up to a subsequence (e n ), (f e (t,x g — v , v)) weakly 
converges to f in L q , where f only depends on (t,x g ,vn, \v±\ — p L ) and is solution of 

^ + «n d x J+ j° Pl e {1 d v J+ (J^Ex) 1 - ■ Vx./ = , (1.10) 

with the initial condition J® (/ )■ 

Moreover, if Jp EE BV{M?), for a.e. p L > 0, then there is no need to extract a subsequence and the 
limit f is the unique solution of (|1.10|) with the initial condition Jp (/ ). 



Remark 1.1 The Bessel operator Jp has some regularization properties. It goes from H s to H s+1 ^ 2 for all 

.1/2 

C X± || /s ±1 

equation (|1.10p . 



s, so that if E G Hi x Hj_ , then Jp E £ H 1 , a condition that ensures the uniqueness of the solution of 



Remark 1.2 The operator J° is important to perform the adaptation of the 6D initial condition to the 5D 
limit model. Indeed, the very fast Larmor gyration creates an initial layer that instantaneously adapts the 
initial condition to the limit model. 

Proof of Theorem 1.1. 

The phase space in position- velocity coordinates is not well adapted to perform the fast gyration averaging. 
To handle it more easily, it is convenient to change the system of coordinates and consider the gyro- coordinates 
defined by 

Xg=X + V ± , Vg = V . (I'll) 

Xg is the position of the gyro-center and p L — \v \ is the ion Larmor radius. To express the gradient in x 
and v in this new system of coordinates, one can conveniently remark that: 

X7 x h(dx g — dVg 1 ) + X7 v hdv g — V Xg hdx g + ~S7 Vg hdv g , (1-12) 

for any smooth function h, with the function h defined by h(x g ,v g ) — h(x,v). Then V x = V Xg and 
V„ = V Vg — V^, where V^ 3 = (d Xg 2 , — d Xg A , 0). Note that V -1 stands for the gradient vector rotated by 
— 7r/2 and not ir/2. The gradients are taken at the corresponding points. For instance, the first equality 
reads \7 x h(x,v) — \7 x Ji(x g ,v g ). The anti-symmetry of _L, a ■ b = —a 1 - ■ b, has been used (and we shall 
make a wide use of it in the sequel). Given these relations, (|1 .4[) can be modified leading one to an equation 
satisfied by the function f s (t, x g ,v g ) — f £ (t, x, v), 



-qJ+ V \\ d *\\ & + E \\ X 9 - V g) <\ fe 



E ± (t,X g - V±) ■ (V v Je - Vif s ) + \v± ■ V v J e = 



(1.13) 



with initial condition / . Here the subscript _L stands for the perpendicular components to the magnetic 
field, for instance Ej_ — (_E l7 i?2,0). Barred quantities are functions of the gyro-coordinates. 



G 



Since (1.6) is a conservative transport equation, i.e. it may be written as d t f + div(. . . ) = 0, the L 9 -norms 
of f e and then of f s are conserved. Thus, ||/e||z,°°(z,£ „) is uniformly bounded. Then, up to the extraction 
of a subsequence, we may assume that f e weakly converges towards some / € L^{L% v ). Upon multiplying 
equation (11.131) by s, in the limit e — > 0, we obtain 

< ■v ,/ = o, 

All the other terms from Eq. (|1.13[) are bounded in the sense of distributions and their product with e therefore 
vanishes in the limit e — > 0. This reduced form of the Vlasov equation implies that the only dependence of 
/ on v 9i ± is on |ug,j_|. Hence with no dependence on the gyrophase. 

Let us now consider equation (j!.13[) for f £ , when integrated against a smooth test fonction <fi with support 
in time avoiding t = (we will handle the initial conditions later) . Let us further assume that the dependence 
on v 9: ± is restricted to a dependence on For such a function, one readily finds that 



feV^Vfidxgdvg = 



for symmetry reasons. As a consequence, the projection of f e takes the following form 

J f e (d t cj) + v gtll d X3Al <f> + E l] d VgAll( l) + E ± - V Vg <P-E ± - V Xg <j>)dx g dv g dt = Q , (1.14) 

keeping in mind that E is calculated at the point x g — v g . We have also used the relation a ■ b = —a • 
as well as the fact that equation (|1.13l) may be written in a conservative form because ( if J is the matrix of 
the linear map v — > v ), 

div Vg (E(t,x-v£)) - Tr((VE)J)--a x2 E 1 +9 Xl E 2 -0, 
div^E- 1 = Tr(J (VE)) = Tr((VE) J) = , 

since E is a gradient. At this stage, we can take the limit and obtain that / also satisfies equation (|1.14l) . 

This is not yet a proper equation in the sense of distributions, even though / depends on Indeed, 
the electric field still exhibits the dependence on the particle position E = E(t,x g — v g ). To obtain an 
equation only depending on |u 9 ,j_|) we use polar coordinates for u s ,j_, i.e. v g — (p L cosip c ,p L sintp c ,vn) and 
Fubini's theorem to integrate the previous integral first in tp Cl then in the other variables. This leads one to 
the following equation: 

f (2irp L )f(d t <l>+v g ,\\d Xgi]l <f> + J° Pl E {] d^+iJ^E) 1 - ■ V Xg <f> 

Jxg,v gA \,r ' ' (1-15) 

+ R(t,Xg,v g ,\\,p L ) ■ Vp L 4>)dXgdv gt \\dp L = , 
where the the operator J" is defined in (jl.8[> . and the term R is equal to 

R(t > x g ,v g<l] ,p L )= f E ± (x g -p L e^ + ^)-e^dp c , 



i) 



with the previously introduced notation e lVc = (cos<p c ,sin<^ c ,0). It can be shown that R is null since it is 
the circulation of the electric field E along a gyrocircle. The latter vanishes in the electrostatic case since 
E is a gradient. With R = 0, (|1.15[) is identical to the equation (jl.101) written distribution wise for 2np L f. 
As there is no dynamics in the p L direction, the factor p L is only a multiplicative constant. One can then 
introduce ^1 p\ +v as test function, and then obtain that fp L /^p 2 L +ri satisfies the equation. Letting i] —> 
in that linear equation, we obtain the appropriate result for the equation governing the evolution of /. 

Regarding the initial conditions we use a similar projection technique with a test function $ depending 
on |u ai|| |, but not vanishing at t = 0. We then obtain equation (|1.14[) . with the right-hand side replaced by 

/ s $(0, x g , v g> n , p L ) dxgdv g: \\p L dp L dLp c . (1.16) 
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The asymptotic limit then leads one to: 

Jh{l f°( x 9,PLe l6 +n p 4) d ^<S>(0,-)dx g dv g ^2Trp L dp L . (1.17) 

Changing coordinate to recover the position-velocity coordinates in the ip c integral, we exactly obtain the 
initial condition J p f° that is expected for /. ■ 



This proof is valid for external electric and magnetic fields. The case where the equation (|1.13[) . invariant 
in the parallel direction, is coupled to the Poisson equation given a Debye length of the same order as the 
Larmor radius has been treated by Frenod and Sonnendriicker in [FS01 . It corresponds to a Debye length 
of order \/e. Technically, we may handle that case with our technic to obtain weak solution as in Arsenev's 
work |Ars75j . The main point is to obtain L p estimates on the density rii = J f dv. They can be obtained by 
classical estimates using upper bounds for the kinetic energy. The case where the Debye length is taken much 
smaller than the Larmor radius will be of greater interest, but the difficult problem is there to average the 
strong oscillations appearing at the scale of the Debye length. We refer to |Gre95j , |Gre96j and |CG00] for 
more details on quasi-neutral plasma without magnetic field, and to [GSR03 for results in the guiding-center 
approximation with p L cx e, Xd oc e and e — > 0. 



1.2 The electroneutrality equation. Heuristic approach. 

In this section we address the self-consistent problem when linking the electric field to the charge distribution. 
The relevant equation is the Maxwell-Gauss equation that relates the divergence of the electric field to to the 
local charge governed by the particle density of charged particles. The latter must be determined using the 
ion distribution function for the guiding centers that is solution of the gyrokinetic Vlasov equation. A similar 
treatment for the electrons must be done. We will follow heuristic arguments together with assumptions 
that are not justified rigorously. However, this heuristic derivation bares some interest. It is an alternative 
to the physicists' presentation of that equation based on the Fourier transform. Furthermore, it allows one 
to recover the electroneutrality equation (|1.25l) used in the GYSELA code to close the gyrokinetic equation 

CLUl. 

For quasineutral plasmas, the electric field response to any charge separation governs a restoring force. 
Should the charge separation extend on a scale larger than the Debye length, the restoring force would be 
to strong to allow any significant charge build-up. The plasma can thus be considered to be everywhere 
with near zero charge, hence quasineutral. As a consequence, the density of negative charges en e (n e being 
the density of electrons) equals the density of positive charges erii (n» being the density of ions). In the 
electrostatic limit, neglecting the time dependence on the vector potential, one can recover this physics based 
argument as an asymptotic limit of the Poisson equation for the electric potential. 

T e n e 

The right hand side is dimensionless and so is e$/T e . The dimensionless control parameter on the left hand 
side operator thus appears as the square of the ratio of the Debye scale divided by the characteristic scale 
of the charge separation. In the limit \ 2 D — > 0, one readily recovers the quasineutrality equation, namely: 
n e = rii. Let us consider the ions density. Given the Vlasov equation and its dependence on the electric 
field, it is likely that the ion distribution function (and thus the ion density) is an implicit function of the 
electrostatic potential. The same applies for the electron density. Then, in the limit X 2 D 0, there is no 
analytical dependence on the electric potential. The latter then becomes a Lagrangian multiplier associated 
to the quasineutrality equation. On such issues, we acknowledge the work of E. Grenier Gre95 , Gre96] (the 
only work, at our knowledge, on that subject where the limiting model is kinetic and not fluid), Y. Brenier 
|BG94j . S. Cordier |CG00j and N. Masmoudi [MasOlj . 

However, if we assume that the system is close to equilibrium, i.e. that the departure from a constant 
electric potential is small, then the electroneutrality equation n e — rii provides an explicit dependence on 
a mean-field electric potential. For an electron-ion plasma, the mass ratio is such that the electrons are 
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far more mobile than the ions. One can then assume that their response to an electrostatic perturbation is 
adiabatic on magnetic field lines or surfaces. In that case, assuming that the equilibrium density of electrons 
n e fi and that of ions n ij0 are equal to no, one may write 

<=_ (#-<#>) e 
n e = n e c «„„(!+ ($_($))), (1.18) 

J e 

where ($) is the average of $ on a closed magnetic field line or surface. The expansion on the right hand 
side holds if e$ << T e . In the very simple geometry that is considered here (with B = (0, 0, B)), (3>) is the 
average in the x\\ direction: ($) = J &(x) dx\\ . The approximation of adiabatic electrons thus reintroduces 
an explicit dependence on the electrostatic potential. 

For the perturbation of the ion distribution, the first difficulty is to obtain the distribution of ions in 
physical space from / written in gyro-coordinates. If the distribution in physical space is assumed to be 
constant on gyrocircles as shown in the last section for the limit model, the following formula is obtained, 

f(t,x,v) = f(t,x + v ± ,\v ± \,v ll ). (1.19) 

Taking the integral in v leads to 

m(t, x)=n J J° Pl (f(t, x, p L , U|| ))2np L dp L dv\\ , (1.20) 

where J° only acts on the x variable (since p L = \v±\ with our conventions). The quasineutrality equation, 
n e = rii, then leads one to 

1 + (*-<*» = J J° PL (f(t,x,p L ,v\{))2np L dp L dv\\. (1.21) 

However, the assumption that the distribution / is constant on the gyrocircles is unrealistic whenever the 
electric potential is not constant. In order to express the inhomogeneity of the density on gyrocircles, we add 
a perturbation to /. It can be expressed as an adiabatic perturbation on the gyrocircles 

-^-^)n Q f t (v), (1.22) 

4> is the average of $ over the gyrocircle of a given particle, and fi(v) is the equilibrium distribution of ions. 
Note that this adiabatic perturbation depends on v± through the choice of the gyrocircle used in determining 
the average $. 

®{t,x,v±) = ^-/$(t,2; + ?r L +p L e i ^)dip c , 
n J $(t, x, v ± )fi(v) dv = ^ J* $(t, x - p L {e l < + e^"))fi{p L , W|| ) dip c dip' c p L dp L dv\\ 

= noJiJVJmx^dp^dp,, (1.23) 

where e ttpc = (cos ip c , sin ip c , 0) and hi(p L ) — 2irp L J fi(p L ,v\\) dv\\. When the equilibrium distribution is a 

maxwellian, fi(v) = 1/ (Tin) 3 / 2 e~\ v \ 2 / Ti , so that hi{p L ) — 2p L /Ti e~ p t^ T \ Adding the adiabatic perturbation 
of the gyrocircles, we finally obtain the following electroneutrality equation, 



l + ± ($_($)) = J Jl L {jtt,x,p„v\\))2irp L dp L dv\\ 

~J (l-{J Q PL f)^x)h t {p L )dp L . 

Multiplying by T e /e, then yields 

($-<$))+! J (<i>-(j;j 2 $) hd fiL )dp L = 



(1.24) 



T 

e 



it. ^ (L25) 
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_ 2 Irp 

Remark 1.3 If the equilibrium density of ions is a Maxwellian, then hi(p L ) = 2p L /Tie l' ' and the 
operator 

r (^ L ) 2 hi{p L )d Pl (1.26) 
is the convolution in the perpendicular plane with the radial function H T (r) defined by 

r 2 

e~ 

H ^ = ^mr- (1 - 27) 

See Appendix\^\ for a detailed proof. 

Since p L is a parameter in the equation of motion (|1.10[) . the equilibrium distribution and the initial 
perturbation remain concentrated on a unique value of p L at all time provided it is the case for the initial 
conditions. A first step to solve the system might be to start under this assumption of single value of p L . This 
would likely be the most difficult step, since the general case is a superposition of such cases. The gyrokinetic 
Vlasov equation (|1.10[) is unchanged when considering a single value of p L , however the electroneutrality 
equation (11.251) may be simplified and becomes: 

($ _ ($) + |(i _ (J^f^( t ,x) = ^(j° Ph - l) , (1.28) 

where rij(f, a;) = 2irp L J fit, x, ) dv\\ is the density of ions in gyro-coordinates. 

Still, after this further simplification, the gyrokinetic Vlasov equation (jl.lOj) coupled to the electroneutral- 
ity equation (|1.28[) remains is a very difficult mathematical problem. The main difficulties are the following. 

• The lack of regularity in the parallel direction. The potential $ has the regularity of / in the parallel 
direction. This is to be compared to the Vlasov- Poisson case where D 2 Q (A$ in the Poisson equation 
given above) has the regularity of /. We may overcome this problem by adding some viscosity in that 
direction, in other words by adding a term —Xd 2 f in equation (|1.10[) . 

• A less important lack of regularity lies in the perpendicular direction. In fact, only the gyro- average 
of <3? appears in (jl.lQI) . Moreover, the term ($) is more regular than <!>, so that by f| 1 . 28[) $ has the 
regularity of (p). Hence J° ($) has the regularity of (J° ) 2 {p)- The operator (J° ) 2 sends L 2 into 
H 1 . That is "almost" enough, in the sense that we could reach a sufficient regularity provided (Jp L ) 2 
were compact from L 2 in H . 

In view of these difficulties, we focus in the next section on the lack of regularity in the parallel direction, 
and consider time-independent solutions depending only on x\\ and v, of the form f(x\u vii)/_l(|i>_l|)- 



2 Steady state solutions in the direction parallel to the magnetic 
field 

Let us consider steady state solutions (/, $) to (|1.10j) - (|1.25|) and let us assume that the function / can be 
written in the form f{x\\ , v\\ )/_l(|«jJ), with J + °° /_l(|w_l|)27t|wj^| d\v±_\ — 1. Then the term /, has no incidence 
in the evolution equation and can thus be factorized in (jl.lQI) . Taking T e = 1 for the sake of simplicity, taking 
again into account the equilibrium density no, not always constant in that section and therefore denoted by 
n, and replacing rii by p, the electroneutrality equation (|1.25[) then reads: 

$-<$>=-- 1. 

n 

where p(z) — J f(z,v)dv. To further simplify the notations, let z = xu and v = v». Solving (ll.10p ~ p.25p is 
then equivalent to finding a distribution function / solution of the following equation. 

df f P\' df 
dz \nJ dv 
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Note that g' is the derivative of g with respect to z. Furthermore, we consider the problem in a slab geometry 
z G [—1, 1], and therefore require given f± as boundary conditions, 

f(-l,v) = /_(»), v > 0, f(l,v) = /+(«), « < 0. (2.2) 



Lemma 2.1 Given n positive, then any solution f to (|2.1|) - (|2.2|l . smc/i </iai £ is non decreasing, must satisfy: 



Proof of Lemma 2.1 The characteristic (Z, V) of (|2.1|) starting from (z,v) is defined by: 

Z'{s) = V(s), Z(0) = z, 

n*)=-(£)W)), no) = «. 

Hence 

l/ 2 ( S )+2^(Z(s)) = v 2 + 2^(z). 
n n 

For u > 0, it crosses {(—1, u),u > 0}. For v < 0, it crosses {(1, it), u < 0} if and only if there is a solution V 
to V 2 + 2£(1) = i> 2 + 2£(ar), i.e. t; 2 > 2^(1) - £(*)). Consequently, 



where 



Hence 



/(*,«) = /_(7(a_ (*,«))) */ v>-J2(?(l)-P(z)), 

\ \n n J 

f(z,v) = f + (V(s + (z,v))) if v< -^ 2 (£(l)-£(zj), 



V 2 {s ± (z,vj) + 2t-(±l) = v 2 + 2^-(z), V(s-(z,v)) > 0, V(s+(z, v)) < 0. 
n n 



/(*,«) = / + (-^2(£ W -£(i)) + ^) i/^-^i)-^) 



Consequently, 



— oo 
+oo 



f-i)l*&)-k-Q)+*y 



V 2 (w( 1 )-sW) 



Changes of variables in both integrals lead to 
For a constant density n, trivial solutions to (I2.1I) - (I2.2I) are 



f(z,v) = f-(v), t»>0, /(«,«) = /+(«), v<0. 
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Proving the existence of solutions to (I2.1[) - (|2.2I) satisfying > 0, for a non constant u, is the aim of this 

section. Denote by Hi, 1 < i < 4, the following set of assumptions. 



(H 2 ) sup / 

z>0 J v 



(H^ ri < 0. 

- - uf (u) 



du =: A < +oo. 



>oJ^ \Ju 2 - x 
(H 3 ) /-(«) = 0, < u < 2 J -^r, where /i = / /+(w)cfo + A. 

V "(1) J-oo 

f-(u)du< f f + (u)du, [ i±^ldu<^. 



Notice that H 2 is satisfied when 



/•-too 

lim uf_(u) < +oo and / it | fL(u) | rfu < +oo. 



Remark 2.4 XTie assumption (H2) ensures that p is a perturbation of the equilibrium n, so that with as- 
sumption (HI) the force-field will always be oriented rightward (no possibilities of trapped particles). 
The assumption (i?3) ensures that there are two beams of ions (one coming from the right and the other from 
the left), and that all the ions coming from one side reach the other side (no turn-back) . 
The (HA) assumption is more technical and enforces that the derivative of the density is bounded. 

Lemma 2.2 Assume Hi, 1 < i < 4, and ^ G L°° . Denote by 

K = {a€ W^d-l, 1]); a > 0, a(l) - a(~l) < Jfc, < a' < 4 M || £ |U} . (2.4) 

There is a solution a £ K to 



»«(*) = f° - 1 " 1 MU) in I r U =LM in 

J-oo yju 2 +a(l)-a(z) JlJl^j V u2 ~ a(z) + a(-l) 



2 



z e [-1,1] . (2.5) 

Moreover, a is the unique non- decreasing solution of (|2.5|) such that 

a(l)-a(-l)G[Q,i]. 

Proof of Lemma 2.2. Prove that the map F that maps a G K in (3 defined by 

n p(z) = f° \u\f + (u) du 
J-oo \/ii 2 + a(l) — a(z) 



2 

+ °° «/-(«) 



^Ju 2 - a(z) + a(-l) 



du, z G [—1,1], 



has a fixed point. First, F maps K in K . Indeed, (3 G M /1 '°°([— 1, 1]) like a, 
f3 is nonnegative, and 



2 



J-oo •/•^a(i)-a(-i) yu 2 — a(z) + a( — 1) 

< T / + ( U )d u + /"tT «/-(u) fe < M , 

J-oo Jy/a{l)-a(-l) \Ju 2 — a(l) + a(— 1) 
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Hence, 0(1) - < ^y. Moreover, 



2^X + ^Y, 
n z n 



where 



x=r . ' » ' - , r - 

J-oo + a(l) - a(z) ^ 2 \J^ri) V u2 a ( z ) + 

y= /° MAW / + °° AW Tdu . 

J-oo (u 2 + a(l) - a(z))i ' 2 \p3b ( m2 ~ + 

Since f3' > 0, and, X < f_ f + (u)du + A = //, then, 

/■+0O /-i \ /■+00 

SU P / / o ,s. f-{ u ) du < -7-nT \ f~(u)du. 



Indeed, either x > -4tt and then 

' - n(l) 



„ _ u 2 n(l) 



(u 2 — x)a 3"\/3£ Sv^A*' 



or x < and then f-(u) ^ implies that u > 2w^j, hence 



n(l) n(l) 



sup / jf-(u)du < — — — g < 

x>0J2^ (u 2 ~x)2 3^3 J_ f + (u)du 



So that, 



and by H4 



Consequently, /?' < Afi \\ ^ W^. And so, F maps K in if. Moreover, F is continuous for the topology of 
G([— 1, 1]), by definition of /3 in terms of a. Finally, F is compact for the topology of G([— 1, 1]), by the 
compact embedding of W 1,0 °([— 1, 1]) in G([— 1, 1]) and the boundedness of K in W 1,00 ([— 1, 1]). It follows 
from a Schauder fixed point theorem that there is a fixed point for F in if. 

Moreover, the solution of (2.5) is unique in the class of non-decreasing functions a such that a(l) — a(— 1) G 
[0, ^j]- Indeed, for any solution a of (|2.5j) in this class, x := a(l) — a(— 1) solves G(x) = 0, where 

G(I) : _ » f" <£M*, + _A_ /° M^W^ 
n(l) JzJ^fc Vu 2 - x n(-l) J_ oa yju 2 + x 

2 /"+oo 



f + (u)du + I f-(u)du 



(i)y-oo »(-i)y 

Using estimates very similar to ()2.6)) and (#4), we can show the function G is increasing. Moreover it 

_2/i_N 



follows from H 3 and i? 4 that G(0) < and G(-$fc) > 0. Hence the value of a(l) - a(-l) is unique, as well 



as (a(l), a(-l)), given by 



«(!) = 4rf /° /+(«)** + ; = = d'U 



a(l) 

W-oo V" 2 + <*(!) -a(-l) io 
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If a and j3 are two non-decreasing solutions of (I2.5[) . then 

{a-0)(z)T(z) = O, ze[-l,l], 



where 



T(z) = ^ 



2 



J -°° v / I1 !-a( J )ttt(l) v /« i -«2)+a(l)(v'« 1 -'>(2)+»(l) + v / » ! ^W+''(l)) 



_ r+oo uj-(u) du 

2 \hBj v u2 - a (' 2 )+ a (- i )\/ n2 -' 3 ( 2 )+ Q (- i )(V M2 - Q ( z )+ a (- i ) + \/ n2 -' 3 ( z )+ a (- 1 )) 

By Si, i? 3 and £T 4 , 

r° /+(«), / +0 ° 



2T(z) > n(l) - / !±^>du - / adu 



>n(l)- f l+^-du-^L [ +0 ° f_(u)du > 0, ze[-l,ll. 

Hence a = /3. 

Theorem 2.3 Assume H t , 1 < z < 4, and f± e L°° . 

There is a unique solution f G £°°([-l, 1] x M) to (|2~T|) - ((2~2|) . such that £ e K. 

Proof of Theorem 2.3 Let a £ K be the solution to (|2.5I) . The distribution function / defined by 

f(z,v) = f-{^/a(z)^a(-l)+v 2 ), v > - ^a(l) - a(z), 
f(z, v) = f+{-yJa{z) - a(l) + v 2 ), v < - a{\) ~ a(z), 

is the unique solution to ([27t]) -(f2~2 |) in L°°([-l, 1] x St) such that % e K. 

A Appendix: The polarization operator as 
a convolution 

In this section we restate and prove the result announced in Remark 11.31 
We define 

F T = ^\ pe-? 2 / T (J° p fdp, (A.6) 



where T is the temperature (of the ions) and J° is the gyro-average operator defined in (|1.8|) . We forget the 
subscript L in the Larmor radius for conveniance. We use here the measure (2p/T)e~ p l T and not an usual 
Maxwellain, because we start from a 2D Maxwellian and perform an integration over the angular variable in 
polar coordinates. 

The operator Fq appears in the electroneutrality equation (| 1 . 25|) . Here we will prove the following propo- 
sition: 

Proposition A.l The operator Fj, is the convolution with the radial function Ht, defined by 

r 2 
4T 



Proof of the proposition. First, the square of the operator J° is 

1 



p 

2-7T p2tt 



JO 



(j;Y(f)(x g ) = — 2 \ I f(x g + pe^ + pe^)dtp c dtp' c 
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To simplifiy it, we use the equality 

p ( e Up e +e Vc) = 2/9 cos 

which helps defining the polar coordinates (r, <fi) 



2p 



cos 



(^) 



(A.8) 



where e = or 7r depending on the sign of the cosine. This is not exactly a change of variable since a couple 
(r,<p) has exactly two pre-images ((p c ,(p' c ) and (ip' c ,ip c ). But since it is always bi-valued, we can use the 
formula for thechange of variables with a factor 2. The presence of e does not introduce specific difficulties, 
if the intervals of integration are distinguished. The Jacobian of this transformation is 



sin 



<Pc - <P C 



so that 



1 



(■£)(/)(*«) = ^ 



2 

2p /-27T 



' 4' 







/(a: 9 + re* 



d(f>dr 



o 



The next step is to introduce cartesian coordinates. This is a standard transform with Jacobian r, so that 

dy 

lB(2p) n^ryj Ap 2 — r 2 

with y e R 2 , r = |y| and the notation .B(a) for the ball of center and radius a. This is exactly the 



(J° P )V)(* S )= I f(x g + y) — 

JB(2p) TT 2 r 



convolution with the radial function 



K(r) 



1(0, 2p) 



(r) 



n 2 r^ 'Ap 2 — r 2 ' 



where x A denote the characteristic function of A. It can be checked that 



hp(r)2nr dr = 1 . 



The last step is to perform the integration in p. As (Jr) is the convolution with the radial function h p , 



7 t = T I P e p ^ T (J p) 2 dp is the convolution with the radial function 



H T {r) = 



In other words, 



H T (r) = 



2 f+°° 
= ^fJ r/2 



h p {r)pe-P 2 ' T dp. 
p 2 /t 



pe 



dp . 



r/2 \J 4p 2 - r 2 
We perform the change of variable p' = \J Ap 2 — r 2 . Hence 

H r(r) = 7^7^ J o e «■ dp 



e 4-t 
2ir 2 rT 

r 2 



2ir 3 / 2 rVT ' 

It can be checked that Ht has total mass one. Indeed, 

POO -j POO 2 

/ H T {r)2Trrdr = —= \ e^^dr=\ 
Jo vvrTjo 
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